Modelos autorregresivos integrados de media móvil

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf 
import statsmodels.tsa.api as smtsa
import statsmodels.tsa.arima.model as arima_model
from statsmodels.tsa.arima.model import ARIMA
import warnings
warnings.filterwarnings("ignore")

url = 'https://raw.githubusercontent.com/evgomez98/wind_speed/main/wind_dataset.csv'
df = pd.read_csv(url)

from statsmodels.tsa.stattools import adfuller
from numpy import log

Prueba de Dickey-Fuller

Iniciamos realizando la prueba de Dickey-Fuller, donde \(H_0\): La serie no es estacionaria.

result = adfuller(df['WIND'])
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])

ADF Statistic: -8.413013
p-value: 0.000000

Con un \(P-Value\)= 0,00 se rechaza \(H_0\). Es decir, la serie de tiempo de la velocidad del viento es estacionaria.

Autocorrelación y Autocorrelación Parcial

def plotds(xt, nlag=365, fig_size=(12, 10)):
    
    if not isinstance(xt, pd.Series):
         xt = pd.Series(xt)
    plt.figure(figsize=fig_size)
    layout = (2, 2)
    
    ax_xt = plt.subplot2grid(layout, (0, 0), colspan=2)
    ax_acf = plt.subplot2grid(layout, (1, 0))
    ax_pacf = plt.subplot2grid(layout, (1, 1))
    
    xt.plot(ax=ax_xt)
    ax_xt.set_title('Time Series')
    plot_acf(xt, lags=365, ax=ax_acf)
    plot_pacf(xt, lags=30, ax=ax_pacf)
    plt.tight_layout()
    
    return None

plotds(df['WIND'])

El primer diagnóstico sobre la gráfica ACF se trata del patrón sinusoidal que presenta, mostrando autocorrelación y sugiriendo estacionalidad en la serie de tiempo. Por otro lado, la función de autocorrelación parcial muestra cortes significativos en múltiples rezagos lo que sugiere que varios términos de autorregresión podrían ser relevantes para el modelo. Debido a la combinación de cortes en rezagos bajos y altos, es posible que sea necesario consideras realizar pruebas de modelos ARIMA con diferentes valores de p y q, también considerar modelos SARIMA que puedan capturar tanto los patrones de autorregresión como los patrones estacionales en los datos.

División del conjunto de datos

Para el modelamiento, se selecciona un año completo para el testeo.

tau_test = 365
tau_val  = 365

train = df['WIND'][:-(tau_val + tau_test)].copy()
val   = df['WIND'][-(tau_val + tau_test):-tau_test].copy()
test  = df['WIND'][-tau_test:].copy()

print(f"Train: {len(train)}, Validation: {len(val)}, Test: {len(test)}")

Train: 5844, Validation: 365, Test: 365

Construcción del Modelo

Para la selección del mejor orden de los componente del modelo, se automatiza a través de un loop que elige con base en el criterio de Akaike.

from sklearn.metrics import mean_absolute_error

best_mae = np.inf }
best_order = None
best_mdl = None

pq_rng = range(1,5)
d_rng = range(0,3)

for p in pq_rng:
    for d in d_rng:
        for q in pq_rng:
            try:
                tmp_mdl = ARIMA(train, order=(p,d,q)).fit()
                val_pred = tmp_mdl.forecast(steps=len(val))
                tmp_mae = mean_absolute_error(val, val_pred)
                
                if tmp_mae < best_mae:
                    best_mae = tmp_mae
                    best_order = (p,d,q)
                    best_mdl = tmp_mdl
            except:
                continue

  Cell In[8], line 3
    best_mae = np.inf }
                      ^
SyntaxError: unmatched '}'

import pickle
with open('best_arima_model.pkl', 'wb') as f:
    pickle.dump((best_order, best_mdl, best_mae), f)
print('hqic: {:6.5f} | order: {}'.format(best_mae, best_order))

hqic: 4.01136 | order: (1, 0, 2)

from statsmodels.graphics.tsaplots import plot_predict

Modelo ARMA

Entrenamos el modelo, segun el mejor orden, en este caso 1, 0, 2. Es decir, un proceso ARMA:

model = ARIMA(train, order= best_order)
model_fit = model.fit()

model_fit.plot_diagnostics(figsize=(14,10));

A través del correlograma se puede concluir independecia en los residuales.

plt.rcParams.update({'figure.figsize': (10,6)})
fig, ax = plt.subplots();
plot_predict(model_fit, 2, ax=ax);
plt.show();

Mediante la función forecast_accuracyautomatizamos el calculo de las metricas de precisión del modelo:

from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from statsmodels.stats.diagnostic import acorr_ljungbox
from scipy.stats import normaltest

index = 'ARMA (1, 0, 2) train'

train_result = ARIMA(train, order=best_order).fit()
train_pred = train_result.predict(start=train.index[0], end=train.index[-1]) 

mae = mean_absolute_error(train, train_pred)
sse = np.sum((train - train_pred) ** 2)
mape = np.mean(np.abs((train - train_pred) / train)) * 100
msd = mean_squared_error(train, train_pred)
r2 = r2_score(train, train_pred)  

ljung_box = acorr_ljungbox(train - train_pred, lags=[90], return_df=True)
   
normality_test_stat, normality_p_value = normaltest(train - train_pred)

df_acc = pd.DataFrame({'MAE': [mae],
                    'SSE': [sse],
                    'MAPE': [mape],
                    'MSD': [msd],
                    'R2': [r2],
                    'Ljung-Box (p-value)': [ljung_box['lb_pvalue'].iloc[0]],
                    'Normalidad (p-value)': [normality_p_value]},
                    index= [index])
df_acc.head()

	MAE	SSE	MAPE	MSD	R2	Ljung-Box (p-value)	Normalidad (p-value)
ARMA (1, 0, 2) train	3.163859	93215.764275	inf	15.950678	0.349119	2.398775e-13	2.992553e-66

index = 'ARMA (1, 0, 2) val'

val = val.reset_index(drop=True)

val_result = ARIMA(val, order=best_order).fit()
val_pred = val_result.predict(start=val.index[0], end=val.index[-1]) 

mae = mean_absolute_error(val, val_pred)
sse = np.sum((val-val_pred) ** 2)
mape = np.mean(np.abs((val-val_pred) / val)) * 100
msd = mean_squared_error(val, val_pred)
r2 = r2_score(val, val_pred) 

ljung_box = acorr_ljungbox(val-val_pred, lags=[10], return_df=True)

normality_test_stat, normality_p_value = normaltest(val- val_pred)

df_acc = pd.DataFrame({'MAE': [mae],
                    'SSE': [sse],
                    'MAPE': [mape],
                    'MSD': [msd],
                    'R2': [r2],
                    'Ljung-Box (p-value)': [ljung_box['lb_pvalue'].iloc[0]],
                    'Normalidad (p-value)': [normality_p_value]},
                    index= [index])
df_acc.head()

	MAE	SSE	MAPE	MSD	R2	Ljung-Box (p-value)	Normalidad (p-value)
ARMA (1, 0, 2) val	3.246131	5883.128472	36.519955	16.11816	0.304795	0.789485	0.000068

index = 'ARMA (1, 0, 2) test'

test = test.reset_index(drop=True)
test_result =  ARIMA(test, order=best_order).fit()

test_pred = test_result.predict(start=test.index[0], end=test.index[-1]) 

mae = mean_absolute_error(test, test_pred)
sse = np.sum((test - test_pred) ** 2)
mape = np.mean(np.abs((test - test_pred) / test)) * 100
msd = mean_squared_error(test, test_pred)
r2 = r2_score(test,test_pred) 

ljung_box = acorr_ljungbox(test - test_pred, lags=[10], return_df=True)

normality_test_stat, normality_p_value = normaltest(test - test_pred)

df_acc = pd.DataFrame({'MAE': [mae],
                    'SSE': [sse],
                    'MAPE': [mape],
                    'MSD': [msd],
                    'R2': [r2],
                    'Ljung-Box (p-value)': [ljung_box['lb_pvalue'].iloc[0]],
                    'Normalidad (p-value)': [normality_p_value]},
                    index= [index])
df_acc.head()

	MAE	SSE	MAPE	MSD	R2	Ljung-Box (p-value)	Normalidad (p-value)
ARMA (1, 0, 2) test	3.325192	6441.950485	60.291396	17.649179	0.355096	0.756981	0.001823

Para el modelo ARMA(1, 0, 2) en el conjunto de prueba, el MAE de 3.33 indica que las predicciones tienen una desviación promedio de 3.33 nudos respecto a los valores observados, lo cual es aceptable en cuanto a precisión. La SSE de 6441.95 revela un error acumulado moderado, mientras que el MAPE de 60.29% muestra que la precisión es limitada en términos porcentuales. La MSD de 17.65 sugiere que el modelo produce errores significativos en sus predicciones. Un R² de 0.36 refleja que el modelo logra capturar una porción modesta de la variabilidad en la serie temporal. Además, el p-valor de la prueba de Ljung-Box (0.757) sugiere que no hay autocorrelación significativa en los residuos, favoreciendo la aleatoriedad de estos. Sin embargo, el p-valor de la prueba de normalidad (0.0018) indica que los residuos no se distribuyen normalmente, lo cual podría impactar la consistencia de los intervalos de confianza y la validez general del modelo. En conjunto, el modelo ARMA(1, 0, 2) ofrece un rendimiento aceptable, aunque con oportunidades de mejora en precisión y ajuste.

import plotly.graph_objects as go
def plot_model(train, val, test, test_pred, title):
  
    train = np.array(train)
    val = np.array(val)
    test = np.array(test)
    test_pred = np.array(test_pred)

   
    train_index = np.arange(0, len(train))
    val_index = np.arange(len(train), len(train) + len(val))
    test_index = np.arange(len(train) + len(val), len(train) + len(val) + len(test))
    test_pred_index = test_index  

    fig = go.Figure()
    fig.add_trace(go.Scatter(x=train_index, y=train, mode='lines', name='Train', line=dict(color='blue')))
    fig.add_trace(go.Scatter(x=val_index, y=val, mode='lines', name='Validation', line=dict(color='orange')))
    fig.add_trace(go.Scatter(x=test_index, y=test, mode='lines', name='Test', line=dict(color='green')))
    fig.add_trace(go.Scatter(x=test_pred_index, y=test_pred, mode='lines', name='Prediction', line=dict(color='red', dash='dash')))

    fig.update_layout(
        title=title,
        xaxis_title="Observaciones",
        yaxis_title="Valores",
        legend=dict(x=0, y=1, traceorder="normal"),
        plot_bgcolor='rgba(0,0,0,0)',
        width=900, height=600
    )

    fig.show()

plot_model(train, val, test, test_pred, 'ARMA Model 1, 0, 2')

La función arima_rolling nos proporciona un bucle para predecir a partir de medias moviles, de modo que el modelo se adapte a la evolución de los datos:

def arima_rolling(train, val, test, best_order):
    history = list(train) + list(val)
    predictions = []
    for t in range(len(test)):
        model = ARIMA(history, order=best_order)
        model_fit = model.fit()
        output = model_fit.forecast()
        yhat = output[0]  
        predictions.append(yhat)
        obs = test[t] 
        history.append(obs)  
        
        print(f'Predicted={yhat:.3f}, Expected={obs:.3f}')
        
    return predictions

Ahora pronosticamos y comparamos con los datos de prueba:

test_pred  = arima_rolling(train, val, test, best_order)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 test_pred  = arima_rolling(train, val, test, best_order)

NameError: name 'arima_rolling' is not defined

index = ARMA rolling (1, 0, 2) test 

mae = mean_absolute_error(test, test_pred)
sse = np.sum((test - test_pred) ** 2)
mape = np.mean(np.abs((test - test_pred) / test)) * 100
msd = mean_squared_error(test, test_pred)
r2 = r2_score(test,test_pred) 

ljung_box = acorr_ljungbox(test - test_pred, lags=[10], return_df=True)

normality_test_stat, normality_p_value = normaltest(test - test_pred)

df_acc = pd.DataFrame({'MAE': [mae],
                    'SSE': [sse],
                    'MAPE': [mape],
                    'MSD': [msd],
                    'R2': [r2],
                    'Ljung-Box (p-value)': [ljung_box['lb_pvalue'].iloc[0]],
                    'Normalidad (p-value)': [normality_p_value]},
                    index= [index])
df_acc.head()

	MAE	SSE	MAPE	MSD	R2	Ljung-Box (p-value)	Normalidad (p-value)
ARMA (1, 0, 2) test	3.326124	6425.376786	61.322275	17.603772	0.356756	0.725282	0.001373

Para el modelo ARMA(1, 0, 2) en el conjunto de prueba, el valor de MAE de 3.33 sugiere una desviación promedio de aproximadamente 3.33 nudos en las predicciones, mostrando una precisión aceptable en términos absolutos. La SSE, con un valor de 6425.38, denota un error acumulado moderado en el ajuste, mientras que el MAPE de 61.32% indica que el modelo no es muy preciso en términos relativos o porcentuales. La MSD de 17.60 confirma una magnitud significativa en los errores, y el coeficiente de determinación \(R^2\) de 0.36 muestra que el modelo logra capturar parte de la variabilidad en los datos, aunque con margen de mejora. El p-valor de la prueba de Ljung-Box (0.725) sugiere una ausencia de autocorrelación significativa en los residuos, favoreciendo la aleatoriedad en estos. Sin embargo, el p-valor de la prueba de normalidad (0.0014) sugiere que los residuos no siguen una distribución normal, lo cual podría afectar la robustez de los intervalos de confianza. En conjunto, el modelo ofrece un rendimiento razonable, aunque podría mejorarse en términos de precisión relativa y ajuste de los residuos.

plot_model(train, val, test, test_pred, 'ARMA rolling Model 1, 0, 2')

pred_error = np.array(test) - np.array(test_pred)

from pandas.plotting import autocorrelation_plot

fig = plt.figure(figsize=(5.5, 5.5))
autocorrelation_plot(pred_error, color='b')
plt.title('ACF for Prediction Error');

Tal como el Ljung-Box no rechaza la independenicia de los residuales del modelo ARMA(1,2), se puede apreciar en la gráfica este comportamiento.

Modelo SARIMA

import statsmodels.api as sm

best_aic = float('inf')
best_pdq = None
best_PDQ = None
best_mdl = None

for p in range(5):
    for d in range(3):
        for q in range(5):
            try:
                tmp_mdl = sm.tsa.SARIMAX(train_data, order=(p,d,q), seasonal_order=(p,d,q,365)).fit()
                tmp_aic = tmp_mdl.aic
                if tmp_aic < best_aic:
                    best_aic = tmp_aic
                    best_pdq = (p, d, q)
                    best_PDQ = (p, d, q)
                    best_mdl = tmp_mdl
            except Exception as e:
                continue

#if best_pdq and best_PDQ and best_mdl:
    #with open('best_sarima_model.pkl', 'wb') as f:
       #pickle.dump((best_pdq, best_PDQ, best_mdl, best_aic), f)

   #print(f'aic: {best_aic:6.5f} | order: {best_pdq} | seasonal_order: {best_PDQ}')
#else:
    #print("No se encontró ningún modelo adecuado.")

El loop para el modelo SARIMA no se ejecuta por el riesgo de sobrecalentamiento del equipo.